# include <stdio.h>
# include <math.h>
# include <stdlib.h>
# define pi 3.1415
# define DU 0.0174533 //
# define rad 57.29578 //
# define JJp 45  //纵波的入射角临界值
# define Rjp 53  //  
# define N 4
# define M 5

/*--------------------------第二次大作业------------------------------*/

/*----------函数定义区---------*/
double snell(double HUDU1, double v1, double v2);
int getMaxinumLaber(double a[N][M], int k);   //获取主元所在行
void R(void);

void main()
{//--------------p是纵波，s横波
	double Vp1=2.134, Vp2=3.018, Vs1=0.860, Vs2=1.595, Q1=2.110, Q2=2.230;
	double V;
	int E; // 入射角度
	double a1, a2, b1, b2; //a1是反射纵波弧度， a2是反射横波, b1是透射纵波， b2是透射横波
	//App纵波反射系数,Aps反射系数，Bpp纵波透射系数，Bps横波透射系数
	int i, j, k; 
	int laber; //计算过渡量
	double sum; //求和
	double temp;
	double x[N];//存放要求的反射系数
	/*-----jjp是通过asin()算出的入射角的临界值，本题目有两个临界角JJp和JJs------------*/
	double APP[JJp];
	double APS[JJp];
	double BPP[JJp];
	double BPS[JJp];
	double a[N][M];//存放值
	FILE *fp1; //App纵波反射系数
	FILE *fp2; //Aps横波反射系数
	FILE *fp3; //Bpp纵波透射系数
	FILE *fp4; //Bps横波透射系数
	
	if((fp1=fopen("App", "w"))&&(fp2=fopen("Aps", "w"))&&(fp3=fopen("Bpp", "w"))&&(fp4=fopen("Bps", "w")) == NULL) 
	{
		printf("erorr!!!");
		exit(0);
	}
	for(E=0; E<JJp; E++)
	{
		/*---角度弧度转化------*/
		V = E * DU;
		/*--------通过snell定律计算----------*/
		a1 = snell(V, Vp1, Vp1);
		b1 = snell(V, Vp1, Vp2);
		a2 = snell(V, Vp1, Vs1);
		b2 = snell(V, Vp1, Vs2);
		/*------赋值看起来有点傻----------*/
		a[0][0] = sin(a1);
		a[0][1] = -cos(a2);
		a[0][2] = -sin(b1);
		a[0][3] = -cos(b2);
		a[0][4] = -sin(a1);
		
		a[1][0] = cos(a1);
		a[1][1] = sin(a2);
		a[1][2] = cos(b1);
		a[1][3] = -sin(b2);
		a[1][4] = cos(a1);
		
		a[2][0] = sin(2*a1);
		a[2][1] = -(Vp1/Vs1)*cos(2*a2);
		a[2][2] = ((Q2/Q1)*(pow(Vs2,2)/pow(Vs1,2))*(Vp1/Vp2))*sin(2*b1);
		a[2][3] = ((Q2/Q1)*((Vp1*Vs2)/(pow(Vs1,2))))*cos(2*b2);
		a[2][4] = sin(2*a1);
		
		a[3][0] = cos(2*a2);
		a[3][1] = (Vs1/Vp1)*sin(2*a2);
		a[3][2] = -(Q2/Q1)*(Vp2/Vp1)*cos(2*b2);
		a[3][3] = (Q2/Q1)*(Vs2/Vp1)*sin(2*b2);
		a[3][4] = -cos(2*a2);
		/*---------------------------*/
		
		//高斯消元法
		for(k=0; k<N; k++)
		{
			laber = getMaxinumLaber(a, k);	
			//主元交换
			if(laber != k)
			{
				for(i=0; i<M; i++)
				{
					temp = a[k][i];
					a[k][i] = a[laber][i];
					a[laber][i] = temp;
				}
			}
			//消元
			for(i=k+1; i<N; i++)
			{				
				if(a[k][k]==0) {		
					break;
				}
				temp = a[i][k]/a[k][k];
				for(j=k; j<M; j++) {
					a[i][j] = a[k][j]*temp-a[i][j];
				}
			}
		}
		//回代求解x
		for(i=N-1; i>=0; i--)
		{
			sum = 0;
			for(j=i+1; j<N; j++)
			{
				sum += a[i][j]*x[j];
			}
			x[i] = (a[i][M-1]-sum)/a[i][i];
		}
		/*-----赋值--------*/
		APP[E] = x[0];
		APS[E] = x[1];
		BPP[E] = x[2];
		BPS[E] = x[3];
	}
	for(i=0; i<JJp; i++)
	{
		fprintf(fp1, "%lf \n", APP[i]);
		fprintf(fp2, "%lf \n", APS[i]);
		fprintf(fp3, "%lf \n", BPP[i]);
		fprintf(fp4, "%lf \n", BPS[i]);
	}
	R();
	if(fclose(fp1)&&fclose(fp2)&&fclose(fp3)&&fclose(fp4) != 0)	
    { 
		printf("the file can not be close!");
		exit(1);
	}
}

//获取主元所在行
int getMaxinumLaber(double a[N][M], int k)

{
	int i;
	int laber=k;
	double maxinum=0;
	
	for(i=k; i<N; i++)
		if(maxinum<fabs(a[i][k]))
		{
			maxinum = fabs(a[i][k]);
			laber = i;
		}
		return laber;
}

double snell(double HUDU1, double v1, double v2)
{
	return asin((sin(HUDU1)*v2)/v1);
}
void R(void) //JD是纵波入射角
{
	double Vp2=2.438, Vp1=3.048, Vs2=1.692, Vs1=1.244, Q2=2.140, Q1=2.400; //原始数据
	//---计算系数-------
	double A;
	double B;
	double Q;
	double z_a;
	double z_b;
	double z_Q;
	//-----------------
	FILE * fp = NULL; 
	double * R_date = NULL;
	int i;
	double a1, b1, a2, b2;
	double V; //弧度
	/*-----------文件打开------------*/
	if((fp=fopen("R_date", "w")) == NULL)
	{
		printf("erorr!!!");
		exit(0);
	}
	
	R_date = (double *)malloc(Rjp * sizeof(double));
	
	Q = (Q1+Q2)/2;
	z_Q = fabs(Q2-Q1);
	
	for(i=0; i<Rjp; i++)
	{
		V = i*DU;
		a1 = snell(V, Vp1, Vp1);
		b1 = snell(V, Vp1, Vp2);
		a2 = snell(V, Vp1, Vs1);
		b2 = snell(V, Vp1, Vs2);
		
		A = (a1+a2)/2;
		B = (b1+b2)/2;
		z_a = fabs(a2-a1); 
		z_b = fabs(b2-b1);
		
		R_date[i]=(1/2*((z_a/A)+(z_Q/Q)) + ((1/2*(z_a/A))-(4*(pow(B,2)/(pow(A,2)))*(z_b/B))-(2*(pow(B,2)/pow(A,2))*(z_Q/Q)))*pow(sin(V),2) + 1/2*(z_a/A)*((pow(tan(V),2))-pow(sin(V),2)));
	}
	for(i=0; i<Rjp; i++)
	{
		fprintf(fp, "%lf\n", R_date[i]);
	}
	
	free(R_date);
	if(fclose(fp) != 0)	
    { 
		printf("the file can not be close!");
		exit(1);
	}
}